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2 : ABSTRACT 
> 

We have derived exact solutions of the isothermal Lane-Emden equation with 
and without rotation in a cylindrical geometry. The corresponding hydrostatic 
equilibria are relevant to the dynamics of the solar nebula before and during 
the stages of planet and satellite formation. The nonrotating solution for the 
mass density is analytic, nonsingular, monotonically decreasing with radius, and 
it satisfies easily the usual physical boundary conditions at the center. When 
differential rotation is added to the Lane-Emden equation, a new class of exact 
solutions for the mass density appears. We have determined all of these solutions 
analytically as well. Within this class, solutions that are power laws or com- 



>• 1 binations of power laws are not capable of satisfying the associated boundary- 

value problem, but they are nonetheless of profound importance because they 
constitute "baselines" to which the actual solutions approach when the central 
boundary conditions are imposed. Numerical integrations that enforce such phys- 



O ' ical boundary conditions show that the actual radial equilibrium density profiles 



O 



emerge from the center close to the nonrotating solution, but once they cross 
below the corresponding baselines, they cease to be monotonic. The actual so- 
lutions are forced to oscillate permanently about the baseline solutions without 
ever settling onto them because the central boundary conditions strictly prohibit 
the matching of the two types of solutions. 

This oscillatory behavior of the isothermal solutions to the Lane-Emden 
boundary-value problem is entirely generic and extends to polytropic models 
as well. Based on our results, we expect that quasistatically-evolving protoplan- 
etary disks should develop oscillatory radial density profiles in their midplanes 
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during the isothermal phase of their collapse. The peaks in these profiles cor- 
respond to local gravitational potential minima and their radial locations are 
ideal sites for the formation of protoplanets; sites that can be accentuated dur- 
ing infall of more matter from the still-collapsing cloud. Indeed, a straightfor- 
ward application to our solar system using an oscillatory solution derived from a 
differentially-rotating baseline yields a highly accurate match (mean relative er- 
ror pa 4%) between the radial density peaks of the model and the semimajor axes 
of all the major planets and the dwarf planets, provided that the mean density 
profile between 0.8 AU and 11.3 AU falls off with radius as R~ 2 5 (corresponding 
to a mean surface density variation of R~ l b that is consistent with the profile 
determined empirically by Weidenschilling in 1977). We believe then that for the 
first time in over two centuries we have a mathematically rigorous explanation 
of all planetary orbits in our solar system and the physics that is responsible for 
planet formation at radii that are not in the least random or arbitrary. 

Subject headings: planets and satellites: formation — planets and satellites: 
general — solar system: formation — planetary systems: formation — planetary sys- 
tems: protoplanetary disks 

1. Introduction 

After a research effort that spans almost ten years, we have been able to derive exact 
solutions of the nonlinear differential equations that describe the equilibrium structures of 
differentially-rotating, self-gravitating fluids with cylindrical symmetry. Our results are di- 
rectly applicable to quasi-equilibrium configurations that may develop early in the evolution 
of protostellar and protoplanetary disks, and we expect that they will help us understand the 
physical conditions that prevail in such systems long before the onset of accretion processes, 
gas ionization, and nonaxisymmetric evolution. The nonrotating analogues of our equa- 
tions have been legendary in the literature; they are known as the nonlinear Lane-Emden, 
Thomas-Fermi, and Emden-Fowler equations. Introduction of rotation to these equations 
complicates matters considerably and very little analytical work has been carried out to date 
in this most interesting case. 

Our results are also relevant to a long-standing problem in the formation of our solar 
system that has preoccupied professional astronomers and nonprofessionals for over two 
centuries, namely the locations of planetary orbits and the observed "order" in the present 
solar system. Naturally, we too have ended up considering the famous Titius-Bode "law" 
of planetary distances but not in an effort to discover some hidden physical principle or 
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previously unknown "universal" underpinning; on the contrary, we saw this "law" as an 
adverse interpretation to different and more complex patterns that we detected in the actual 
data and that we describe in § 11.31 We begin in §§ I1.1H1.2I with a brief overview of the 
research that has been carried out by other researchers on these topics. We also mention in 
§ ll.3l some of the history of our work that has finally tied up together all the pertinent issues. 
We refer especially to our past unsuccessful attempts to resolve the question of planetary 
distances because it is such setbacks that eventually pave the way for a successful conclusion 
of the research effort. 



1.1. The Lane— Emden Equation 

The Lane-Emden differential equation (Lane 1869-70; Emden 1907) describes the equi- 
libria of nonrotating fluids in which internal pressure balances self-gravity. Spherically sym- 
metric solutions of this equation came to the attention of astrophysicists when Chandrasekhar 
included them in his 1939 monograph "An Introduction to the Study of Stellar Structure," 
but interest in such solutions continued to be largely academic because real stars rotate and 
rotation destroys spherical symmetry and modifies their internal profiles and physical char- 
acteristics. In the latter half of the twentieth century, the isothermal solution, commonly 
referred to as the "singular isothermal sphere," and its nonsingular modifications found some 
interesting applications to the structures of collisionless systems such as globular clusters and 
early-type galaxies (Binney & Tremaine 1987; Rix et al. 1997); to the structures of large- 
scale gaseous systems such as X-ray halos around elliptical galaxies (Fabbiano 1989) and 
clusters of galaxies (Sarazin 1988; Fabian 1991); and to simplified models of gravitational 
lenses (Kochanek 1995). 

Emden's work also attracted the attention of physicists outside the field of astrophysics 
(Fowler 1914; Thomas 1927; Fermi 1927) who studied generalized polytropic forms of the 
Lane-Emden equation for specific polytropic indices n. The extensive studies of Fowler 
(1914, 1930, 1931) produced some singular solutions for n = 3 and established the so-called 
Emden-Fowler equation in the literature, while the works of Thomas (1927) and Fermi (1927) 
produced the so-called Thomas-Fermi equation, an important milestone in atomic theory. 
At present, both of these equations continue to be the subjects of investigations by physicists 
and mathematicians alike. Physicists are drawn to the Emden-Fowler equation because it 
appears in the kinetics of Landau-Ginzburg critical phenomena (see the detailed account 
of Dixon & Tuszyhski 1990) and in the kinetics of combustion (Frank-Kamenetskii 1955). 
Mathematicians use these nonlinear equations as laboratories to study a wide variety of 
properties in their solutions — positivity, uniqueness, singularities, monotonic vs. oscillatory 
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behavior, and bifurcations — usually without having analytic solutions at hand (e.g., Wong 
1975; Lions 1982; see also Goenner & Havas [2000] for a comprehensive list of all known 
exact analytic solutions; as well as Berkovich [1997] and Goenner [2001] for classifications of 
some solutions obtained from Lie-group symmetries). 

The cylindrical form of the classical Lane-Emden equation differs from the spherical 
form in just one coefficient and both forms have been included as special cases in studies of 
the generalized Emden-Fowler equation (Horedt 1986; Dixon & Tuszyhski 1990; Berkovich 
1997; Goenner & Havas 2000). But until now, the cylindrical Lane-Emden equation has 
not attracted attention on its own merit — not even in astrophysics where it can be used 
as a simplified (but nonlinear) model of a disk-like or cylindrical self-gravitating gas. The 
only analytic studies of the cylindrical form that we are aware of have been made by Jeans 
(1914) and Robe (1968) who found a Bessel-function solution of the linear n — 1 case with 
uniform rotation and by Stodolkiewicz (1963) and Ostriker (1964) who solved the nonrotating 
isothermal (n — > oo) case. Ostriker (1964), in particular, also showed that the radius and 
the mass per unit length are finite in all models with < n < oo and that, in contrast to 
the singular isothermal sphere, the nonrotating, nonsingular, isothermal cylinder has finite 
mass per unit length in spite of its infinite radial extent. 

After the work of Robe (1968), some numerical studies of rotating isothermal cylinders 
and disks (Hansen et al. 1976; Schmitz & Ebert 1986; Narita et al. 1990) and rotating 
polytropic cylinders (Schneider & Schmitz 1995) also found nonmonotonic solutions for the 
equilibrium density akin to the Bessel-function oscillatory solution of the n = 1 polytrope in 
uniform rotation; but understanding this behavior proved elusive and so the nonmonotonic 
density profiles were marginalized. Schneider & Schmitz (1995), in particular, adopted some 
differentially rotating profiles and obtained numerical solutions for the radial density profiles 
of models with negative polytropic indices (n < —1). These authors did not study the 
isothermal models of interest here, and they did not explain the unusual properties of their 
solutions. Nevertheless, their work is by far the most closely related study to our work and 
their polytropic results should be considered in conjunction with our results from isothermal 
models. 



1.2. The Titius-Bode Algorithm 

The numerical algorithm called the Titius-Bode "law" has been known for 240 years 
(e.g., Nieto 1972; Lecar 1973; Danby 1988). It relies on an ad-hoc geometric progression to 
describe the positions of the planets in the solar system and works fairly well out to Uranus 
but no farther (Jaki 1972). The same phenomenology has also been applied to the satellites 
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of the gaseous giant planets (Neuhauser & Feitzinger 1986). Two modern brief reviews of the 
history along with criticisms of this "law" have been written by Graner & Dubrulle (1994) 
and by Hayes & Tremaine (1998). Currently, the general consensus is that a satisfactory 
physical basis has not been found for this numerical coincidence despite serious efforts by 
many researchers in over two centuries. Furthermore, opinions differ on whether such a 
physical basis exists at all. 

Apparently, many researchers still believe that the Titius-Bode algorithm does have a 
physical foundation and continue to work on this problem. The last decade of the twentieth 
century, in particular, saw a resurgence of investigations targeting precisely two questions: 
the origin of the "law" (Graner & Dubrulle 1994; Dubrulle & Graner 1994; Li, Zhang, & 
Li 1995; Nottale, Schumacher, & Gay 1997) and its statistical robustness against the null 
hypothesis (Hayes & Tremaine 1998). 

Hayes & Tremaine (1998) found some statistical evidence that the Titius-Bode "law" 
could be related to the long-term dynamical stability of the solar system. Their work is 
the latest in a long line of differing arguments made by several statisticians in the past and 
summarized by these authors (see also Lynch 2003). More importantly, Hayes & Tremaine 
(1998) dismissed without discussion all the previous purported explanations of the physical 
origin of the "law" as "not entirely convincing." In our opinion, that is an example of good 
physical instincts, but it still requires some physical understanding of planetary distances 
before final judgment is passed. We believe that the work presented below does provide the 
required understanding, and we will return to a broad discussion of the Titius-Bode " law" 
in § [472] below. 

1.3. Past Analyses of Planetary Distances and the Lane— Emden Equation 

Just as Jeans (1914) had done years ago, we began studying rotating self-gravitating 
cylinders in 1997 for the same reason; as Jeans put it: 

"All the essential physical features of the natural three-dimensional problem appear to be 
reproduced in the simpler cylindrical problem, so that it seems legitimate to hope that an 
argument by analogy may not lead to entirely erroneous result [s]. " 

Being unaware of the work published in French by Robe (1968), we solved analytically 
the n — 1 polytropic Lane-Emden equation with cylindrical symmetry, uniform rotation, 
and proper boundary conditions (Christodoulou 1997, unpublished). (Applying boundary 
conditions did not concern Jeans who was interested in the stability of local deformations in 
compressible gases.) The radial density profiles for n — 1 and for sufficiently fast rotation 
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Table 1: 

Planetary Distances in Our Solar System 



Index 


Object 


Semimajor 


Titius-Bode 




Inversion 




Mean 
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Name 


Axis Oj 


Distance 


Error 


Distance 


Error 


Distance 


Error 


Mi 






(AU) 


(AU) 


(%) 


(AU) 


(%) 


(AU) 


(%) 




1 


Mercury 


0.387 


0.4 


3.4 












2 


Venus 


0.723 


0.7 


-3.2 


0.622 


-14.0 


0.694 


-4.1 


1 


3 


Earth 


1 


1.0 


0.0 


1.050 


5.0 


1.124 


12.4 


2 


4 


Mars 


1.524 


1.6 


5.0 


1.663 


9.1 


1.883 


23.5 


2 


5 


Ceres 


2.765 


2.8 


1.3 


2.816 


1.8 


3.364 


21.6 


2 


6 


Jupiter 


5.203 


5.2 


-0.1 


5.135 


-1.3 


6.151 


18.2 


2 


7 


Saturn 


9.537 


10.0 


4.9 


9.992 


4.8 


12.20 


27.9 


2 


8 


Uranus 


19.19 


19.6 


2.1 


16.93 


-11.8 


19.80 


3.2 


1 


9 


Neptune 


30.07 


38.8 


29.0 


27.52 


-8.5 


29.34 


-2.4 


1 


10 


Pluto 


39.48 


77.2 


95.5 
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Magnification Ratio: 
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- di-i). 











exhibited permanent oscillations due to the zeroth-order Bessel function Jo(R) that is the 
dominant part of the solution. Numerical integrations for polytropes with n > 1 also showed 
that the solutions for the density oscillate permanently with radius. This generic oscillatory 
behavior then led us to consider the distribution of semimajor axes of planetary orbits in 
the solar system and their possible connection to the radial density peaks found in the 
equilibrium solutions. 

On closer inspection of the actual planetary data, we found patterns other than the 
Titius-Bode algorithm that could in principle also provide good fits to some sections of the 
data. These patterns are summarized in Table 1 where we list the observed semimajor axes 
ctj of planetary orbits in our solar system (e.g. Kaufmann 1994) along with three empirical 
fits to the data, the Titius-Bode "law," inversion, and an arithmetic mean. The inversion 
distance for each orbit is the geometric mean of the actual semimajor axes of the two nearest 
neighbors. The mean distance for each orbit is the arithmetic average of the actual semimajor 
axes of the two nearest neighbors. Relative errors are calculated for all three fits with respect 
to the observed values. The inversion distances are listed in Table 1 because this is the first 
pattern that we saw in the observed data rather than the Titius-Bode "law." The mean 
distances are shown for comparison purposes; they are very accurate only for those orbits in 
which either or both of the other two fits fail. 
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The surprisingly small errors in the arithmetic-mean distances of Venus, Uranus, and 
Neptune led us to the working hypothesis that the effect that may cause inversion in the 
intermediate orbits does not operate at small or at large distances; moreover, it is smoothly 
replaced at the two ends by a slick new regularity that manifests itself as a modest arithmetic 
average. This smooth transition from inversion to arithmetic averaging was another signifi- 
cant pattern that we saw in the data: For example, the orbit of Uranus (along with Jupiter's 
orbit) reproduces successfully the orbit of Saturn by inversion while, on the other side, the 
Saturn-Uranus pair is clearly in arithmetic progression with Neptune; and the Earth- Venus 
pair shows inversion on the side of Mars and arithmetic progression on the side of Mercury. 
We thought that this could not be a numerical coincidence because we found similar smooth 
transitions in the orbital distances of the regular (and even some irregular) satellites of the 
Jovian planets. 

The geometric-mean spacing of the 1-10 AU objects in Table 1 implies that their orbits 
are inverted images of every other one with respect to the corresponding in-between orbit 
(e.g., Coxeter 1989), i.e., that the semimajor axes of any three consecutive orbits obey the 
relation 

Oi_i • a i+1 = a] , (3 < i < 7) . (1) 

In geometry, inversion is the gateway to hyperbolic space, where the inverted orbits would 
appear to be conveniently equidistant — a perfect symmetry indeed, taking place in a space 
that we cannot even visualize. But we could not see how the orbital plane of the solar 
system could be so strongly curved in its middle and so flat at the two ends; so eventually 
we abandoned this line of reasoning and the high degree of symmetry in hyperbolic space. 

Next we turned to geometric optics, the only part of physics where the same relation 
occurs. Eq. ([1]) can also be written as 

1 1 1 

a>i — aj + i — cii <ii 

This is a mirror equation and implies that the (i — 1) and the (i + 1) orbits are mirror 
images of one another while the in-between orbit plays the role of a concave mirror. Also, 
the magnification ratio 

Mi = , (3) 

is clearly larger than 1 and approximately constant throughout the system of orbits in ge- 
ometric progression (i.e., for 3 < i < 7); and it reduces to M.{ ~ 1 for those orbits in 
arithmetic progression (see Table 1). This model appeared promising for a while because it 
suggested that some mechanism could potentially be responsible for tapering off the geomet- 
ric progression to an arithmetic progression, as indicated by the magnifications of the orbits. 



(2) 
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Unfortunately, the mirror equation of geometric optics is not based on a fundamental prin- 
ciple per se, it is only a linear approximation valid for paraxial rays, and we did not think 
that its deeper underlying principle (Fermat's principle) could be applicable to planetary 
orbits. Because a mirror equation similar to eq. (j2j) is not derived from first principles in 
any other part of physical theory, we finally became convinced that the observed planetary 
distances, the Titius-Bode algorithm, and the inversion distances could not be derived from 
a physical principle; and we redirected our effort toward mathematical reasons that might 
be responsible for the observed order in the solar system. This prompted us to return to the 
Lane-Emden equation and to focus exclusively on its intrinsic properties in the cases with 
and without rotation. 

In the case of no rotation or for some specific rotation profiles, the Lane-Emden equation 
is scale invariant and can be transformed to an autonomous differential equation (see e.g. 
Bender & Orzag [1978] and Visser & Yunes [2003]) for the theory and the transformations of 
such differential equations). Scale invariance was also exploited by Graner & Dubrulle (1994) 
and Dubrulle & Graner (1994) who argued that cold, self-gravitating, perfect-fluid disks are 
scale invariant and that this condition is sufficient to generate unstable radial modes that 
are equidistant in lni? (see also Schmitz [1984] and Li et al. [1995] for stability analyses 
that have effectively led to the same result). Unfortunately, the scale invariance of the 
inviscid hydrodynamical equations is easily broken by the chosen boundary conditions for 
the equilibrium system and for the stability problem; because of this fact, we thought that we 
should not look for a geometric progression of the Titius-Bode kind in the stability problem. 
After all, scale-invariant unstable modes, such as those studied by Dubrulle & Graner (1994), 
could only produce a generic geometric progression and they would be incapable of matching 
the observed loss of inversion at small and large distances, where planetary orbits seem 
to taper off neatly to two different arithmetic progressions (see Table 1). Therefore, a 
mechanism based on such modes of disturbance would suffer from the same problem that 
also afflicts inversion to a hyperbolic space, the Titius-Bode geometric progression, and any 
other idea that overemphasizes the geometric spacing of the intermediate planetary orbits 
and ignores the observed arithmetic progressions of orbits in the inner and the outer solar 
system. 

Returning to the scale invariance of the Lane-Emden equation itself, we realized that it 
would not matter if this equation lost this symmetry by the applicable boundary conditions so 
long as its autonomous form were to exhibit discrete scale invariance (DSI; Sornette 1998). 
In the theory of autonomous systems, DSI is a stronger constraint than scale invariance, 
and it is associated with limit cycles in the phase portrait of the differential equation (see 
Appendix B in Visser & Yunes [2003]). A solution that exhibits DSI is not invariant for 
any arbitrary rescaling of the independent variable, but it still is self-similar for a specific 
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rescaling of that variable. This fact led Sornette (1998) to propose that the Titius-Bode 
"law" may be a manifestation of DSI in the solar system, where the constant ratio of the 
Titius-Bode geometric progression or, equivalently, the constant magnification ratio M.i ~ 2 
of inversion is the specific factor that rescales each planetary orbit to the next outward orbit. 
We followed then the hypothesis of Sornette (1998) and the methodology of Visser & Yunes 
(2003) and we constructed numerically the phase portraits of the autonomous forms of the 
polytropic and isothermal Lane-Emden equations, where we looked for limit cycles and 
discrete self-similar behavior with absolutely no success. 

At that point, being convinced that there is no fundamental physics behind the observed 
patterns shown in Table 1 above, and with the theory of nonlinear autonomous systems re- 
vealing no intrinsic periodicities, we basically had only one remaining option: to return 
to our initial approach and solve the full boundary- value problem for the differentially- 
rotating equilibrium systems in order to see if any interesting patterns would emerge from the 
differentially-rotating solutions — patterns similar to the Bessel-function oscillations found 
in the uniformly-rotating n = 1 polytropic case that did not exist intrinsically in the differ- 
ential equation but were generated and governed by the applied boundary conditions. We 
proceeded to do just that, and our results for the isothermal case are described in the sections 
that follow below. In the end, we were truly surprised by the simplicity of the mechanism 
that generates the patterns seen in Table 1; because the discord between the equilibrium 
density profile favored inherently by the differential rotation and the profile imposed exter- 
nally by the boundary conditions, although low-key and inconspicuous, it is nonetheless a 
plain fact and does provide a simple resolution to the long-standing, centuries-old problem 
of planetary "order" in our solar system (see § 12. 4p . 

1.4. Outline 

The remainder of the paper is organized as follows. In § |2l we solve analytically for the 
equilibrium structure of the midplane of a gaseous isothermal disk, incorporating in the Lane- 
Emden equation the effects of self-gravity, differential rotation, and thermal pressure. In § [3j 
we adopt a four-parameter analytic solution as our baseline and we use the rotation profile 
of the baseline to compute numerically the corresponding oscillatory equilibrium solution 
that obeys physical boundary conditions at the center. Then we fit the density maxima of 
this solution to the planetary orbits in the present solar system in order to determine the 
underlying physical characteristics and the stability properties of the baseline model. In § HI 
we discuss the significance of our results for our solar system and for protoplanetary disks 
in general. 



-10- 



2. Isothermal Equilibrium Models 

We consider the axisymmetric equilibria that are available to a rotating self-gravitating 
gas in the absence of viscosity and magnetic fields. We adopt cylindrical coordinates (R, 0, z) 
and the assumption of cylindrical symmetry (d/dz = 0) which lets us ignore ^-dependent 
gradients and reduces the problem to one dimension, the distance R from the rotation axis. 
This technique has become common practice in studies of rotating, self-gravitating, fluid 
disks (e.g., Goodman & Narayan 1988; Christodoulou & Narayan 1992; Christodoulou 1993; 
Christodoulou, Contopoulos, & Kazanas 1996) because it simplifies the stability analyses 
of effectively two-dimensional modes of disturbance. In what follows, we are interested 
in equilibrium structures that describe the physical conditions across the midplane of a 
protoplanetary disk, so the assumptions d/dcf) = d/dz = allow us to tackle the problem by 
solving ordinary differential equations (ODEs). 

We further adopt a rotation law of the form 

Q(R) = Q f(x) , (4) 

where x = R/Ro is a dimensionless radius and the length scale Rq will be specified in eq. (jUJ) 
below. Furthermore, Q(R) is the angular velocity, Qq is the value of Q at some fixed radius, 
and the dimensionless function f(x) for differential rotation is generally an arbitrary function 
of x. For centrally condensed models, it is convenient to choose f2 = O(0) and the regularity 
condition /(0) = 1, while for singular or annular models, we choose Qq = Q(Rq) and the 
normalization /(l) = 1. 

Finally, we assume an isothermal equation of state of the form 

P = clp , (5) 

where P is the thermal pressure, p is the gas density, and Cq is the constant isothermal 
sound speed. Additional calculations in which we used polytropes with indices n = 1 — 3 
and the results of Schneider & Schmitz (1995) who used polytropes with indices n < — 1 
demonstrate that the characteristics of the solutions presented below are largely insensitive 
to the particular choice of n. 



2.1. The Lane— Emden Equation With Rotation 

Axisymmetric and cylindrically symmetric, nonmagnetic equilibria for a perfect fluid 
are described by the equation 
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where the gravitational potential <&(R) satisfies Poisson's equation 

where G is the gravitational constant. Combining eqs. (j4])-(17]) and using the definition 
x = R/Ro, we find a second-order nonlinear innonhomogeneous ODE that can be cast in 
the form 

1 ^#lnr + r = ^(x 2 f) , (8) 
x ax ax 2x ax 

where r = p/po, po is the maximum density or a fixed cutoff density for singular models, 

(3 = Qq/Qj, Q 2 j = 2irGpo, and 

r 2 r 2 



The term f2j represents the gravitational (Jeans) frequency and the dimensionless rotation 
parameter p measures centrifugal support against self-gravity; in general, < Po < 1, since 
the gas is also partially supported by pressure gradients in the radial direction. 

Eq. ([S]) reduces to the classical isothermal Lane-Emden equation in the absence of 
rotation (P = 0)0 We derive analytically the nonrotating solutions in § 12.21 using the 
modern theory of nonlinear ODEs (e.g., Bender & Orzag [1978]; see also the classic works 
of Stodolkiewicz [1963] and Ostriker [1964]). Then, in § 12.31 we derive analytically a class 
of particular solutions of the full problem (eq. [8] with arbitrary f(x) differential rotation) 
and, in § \2A\ we discuss a subset of composite power-law solutions that are astrophysically 
interesting despite the fact that they are incapable of obeying the proper boundary conditions 
at x = 0. 



2.2. Nonrotating Solutions 

In the absence of rotation, the isothermal Lane-Emden equation (eq. [8] with p = 0) 
reads 

1 d d , , 
-— x— lnr + r = 0. 10 

x ax ax 

This equation is scale invariant under the transformation [x — > A x, r — > r), where A is 
an arbitrary constant and p = —2. Therefore, it can be transformed to an autonomous form: 



1 Ot for a flat rotation curve with f(x) — 1/x. This singular rotation law is not of interest in this work. 
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Using this value of p, we define r = x~ 2 w(x) and we cast eq. ffTU]) into an equidimensional- 
in-a; equation for the new function w(x). This equation is invariant under the transformation 
x ^ X x (i.e., its p — 0) and reads 

d 2 d 

x 2 — -law + x— law + w = . (11) 

dx dx 

Finally, we let y(x) = law and we Euler-transform the independent variable [x = e ) to 
obtain the autonomous form for the function y(t): 

y + eV = 0, (12) 

where the dots denote derivatives with respect to t. Since the first derivative is missing from 
eq. ( EEZD , we can integrate twice. The first integral is 



f = Cx _ 2e y } (is) 

and the solution is 

C 2 ±t, (14) 



where C\ and C 2 are integration constants. The integral in eq. (|14p can be calculated in 
closed form 

<uj _ i :ln v^Zg7g-i ; (15) 



VCi - 2eV VC~i v/1 - 2 e y/d + 1 
and a series of backsubstitutions produces the following general solution for the density r{x): 

t(x) = 2Ak 2 , X - j-r , (16) 
v ; (1 + Ax k ) 2 v ; 

where A and k are arbitrary positive constants. The condition that A > is physical 
and ensures that the density profiles are nonnegative. The condition that k > is not a 
limitation: eq. ffl6|) is invariant under the transformation (k — > —k, A — ► 1/A) because k 
contains implicitly the ± duality seen in eq. fTl4l above. Therefore, only positive values of k 
need to be considered without this causing loss of generality. 

The equilibrium solutions obtained in eq. ( TT6T) can be classified into three types: 



1. For < k < 2, the density profiles are singular at x = and decrease monotonically 
for x > 0. These solutions are analogues of the singular isothermal sphere. 
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2. For k = 2, the density profile is centrally condensed and can easily satisfy proper 
boundary conditions at x — 0. Using A = 1/8, eq. ffTB"]) reduces to the Stodolkiewicz- 
Ostriker solution 



(l+X 2 /8) 2 



for which r(0) = 1 and ^(0) = 



3. For k > 2, each density profile has a hole at x = — i.e., r(0) = — and peaks at a 
finite radius x*, namely 

= (\-^) 1/k ■ (18) 



.A £; + 2, 

where r(x*) = 1 and ^(^*) = 0; this normalization fixes the value of A for any choice 
of k > 2. The two constants are related by the equation 



A 2 {k-2) k ~ 2 {k + 2) k+2 = 2 k . (19) 

The mass in all of these models is strongly concentrated around the density maximum 
x = x*; this gives them the appearance of slender annuli despite the presence of 
extended regions of very low densities on either side of the maximum; regions that 
extend all the way to x = and to x = oo. Also note that the solutions for k > 3 
are concave and shallow near the center where they have, not only r(0) = 0, but also 
^(0) = 0; while the solutions for 2 < k < 3 are convex and steeply rising near the 
center where they have dr/dx — > 2Ak 2 (k — 2)/x 3 ~ k as x — > 0. 

At large radii (x » 1), all solutions are decreasing with radius and the density falls off as 
x -k-2 Thjg rapid asymptotic decline, which is steeper than x~ 2 , is responsible for keeping 
the mass per unit length /i finite in all models. A direct integration of eq. fll6p shows that 
/i/27rp -Ro = independent of the constant A. Letting k = 2 in this equation, we recover 
Ostriker's (1964) result, /i = 8np Rl, for the centrally condensed cylinder. 



2.3. Rotating Solutions 

The isothermal Lane-Emden equation with differential rotation (eq. [8]) is repeated here 
for the purpose of discussion: 

~x^r + r = ^(x 2 f 2 ) . (20) 

When the right-hand side (hereafter RHS) of this equation is nonzero (i.e., when /3 ^ 
and f(x) 7^ 1/x), the property of scale invariance is lost from all cases of interest (uniform 
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rotation, power-law rotation, etc.), irrespective of the prescription chosen for the differential 
rotation taction /(x)H 

Eq. (1201) has no special symmetry associated with it, and this is 
enough to make many researchers turn the other way. This is probably why some interesting 
features of eq. f[2"0"j) that we describe below have gone unnoticed for so long. 

The RHS of eq. (1201) is not merely a rotation-dependent correction term to the classical 
isothermal Lane-Emden equation (flQj) . The introduction of rotation changes the properties 
of the ODE to such a large extent that the nonrotating solutions found in § 12.21 cannot guide 
the effort to find rotating equilibrium solutions. In fact, it is the functional form of the 
RHS that determines now the structure of the solutions of the entire ODE: By inspection of 
eq. (1201) . we can write down an entire class of particular equilibrium solutions, namely 

= is (* 2/2 > • (21) 



provided that 



d d , . . 

— x— lnr = , (22) 
dx dx 



also holds true. Using eq. (l2Tj) . we write 

lnr = ln^f - lnx + ln-^-(x 2 / 2 ) , (23) 
2 dx v ' 

and substituting this form into eq. (1221) we find an ODE for all the differential-rotation laws 
f(x) that satisfy eq. (J2"2"j) identically and make eq. (j2"Tp a family of exact solutions of the 
Lane-Emden equation with rotation (eq. [20]): 

-^-fln-f (* 2 / 2 ) = 0. (24) 
dx dx dx 

This third-order linear ODE can be integrated directly to yield the following results: 

j- ^f) = Ax k , (25) 



implying that 



r{x) = 4- ■ Ax k ~ x , (26) 



2 Eq. PD)) with a nonzero RHS is scale invariant only for f(x) — \J A In x + B/x, where A and B are 
arbitrary constants. This case can be solved by transforming the scale-invariant ODE to its autonomous 
form (as was explained in § [272]) , but there is no need to do so; its solution is obtained easier by the method 
described in this subsection. 



- 15 - 



and 



/(*) 



y/Ag(x) + B 



(27) 



x 



where A, B, and k are arbitrary integration constants (that are unrelated to those used in 



implying that dg/dx = x k for all values of k. With so many free parameters (A, B, and k) 
in the differential-rotation profile, these results can easily become a theorist's playground. 
Here we highlight just a few interesting points: 

1. Parameter Constraints. — Eq. ( 1261) shows that r(x) > only for A > 0. This constraint 
also limits the physical values of k when B < in eq. ()27[) ; for example, k > — 1 when 
B = 0. This limitation can be easily circumvented by implementing composite rotation 
profiles with B > (see item 4 in this list and § 12.41 below). 

2. Monotonically Decreasing Profiles. — Eq. f[2"6"j) shows that r(x) is a decreasing function 
of x for | A; | < 1. The same condition is sufficient to also make f(x) a decreasing 
function of x provided that B > in eq. (!27j) . 

3. Uniform Rotation. — For A = 2, B = 0, and fe = 1, eq. (1271) reduces to /(x) = 1 and the 
equilibrium density (eq. [26]) then is r(x) = (3q = constant. Note that this constant 
cannot be adjusted freely, e.g., it cannot be reset to 1; the requirement that f(x) = 1 
fixes A in eq. ( 1271) and then the uniform density gets fixed to /3q by eq. (1261) . 

4. Composite Profiles With B > 0. — Steep density profiles with k < — 1 can be obtained 
by selecting B > and by incorporating a central region of uniform rotation (see § 12.41 
for details). Even more complex equilibrium profiles can be constructed by combining 
two or more power laws. A composite profile is demonstrated in § 13.11 below, where 
we connect two disjoint regions of constant density with a power law and we apply the 
result to the early structure of the solar nebula. 

5. Asymptotic Regime. — For k < — 1 and B > 0, eq. (T28|) shows that g(x) — > as x — » oo 
and eq. (T27|) then exhibits the asymptotic behavior f(x) — > \fBjx. Therefore, all steep 
density profiles with k < — 1 and r(x) oc x k ~ 1 approach a flat rotation curve (QR — > 
constant) at large radii, independent of the value of k. 



I2~2l) and 




(28) 



From the perspective of the physics that dictates the above profiles, the solutions (I2T1) 
of the Lane-Emden equation (!20l) describe a class of rotating self-gravitating equilibria in 
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which z-gradients are neglected and the radial gradient of the gravitational acceleration is 
balanced exactly by the radial gradient of the centrifugal acceleration at every radius R. This 
occurs because, in the Lane-Emden equation, we have gone to second order by taking an 
extra derivative on the components of the equation of hydrostatic equilibrium. The balance 
of gradients can be seen, most easily, by substituting eq. (12T|) into the one-dimensional 
Poisson's equation V 2 ?/> = r, where if) = $/cq is the normalized potential; the result is 



1 d 

x dx 



x 



dip 
dx 



1 d 

x dx 



(29) 



In this equation, the bracketed terms are the gravitational and centrifugal accelerations, 
respectively. This type of balance is different than the balance commonly discussed between 
the magnitudes of these two accelerations in rotating gravitating systems; and the power-law 
density solutions are borne out of this conformance of the two gradients, whereas the familiar 
stalemate between centrifugal and gravitational force is only relevant to purely homogeneous 
fluid equilibrium systems or particle systems with no pressure support. In the isothermal 
gaseous case of interest here, a power-law density profile satisfies naturally the condition 
that the radial variation of the enthalpy gradient p~ l dP/d\nR be zero (see eq. |22j ) and 
so the pure power-law profile is not at all influenced by the radial variation of the pressure 
gradient — it is an exact solution of eq. ff20l PI 



2.4. Composite Models and Boundary Conditions 

Many of the rotation profiles discussed above are singular at x = 0. The solutions for 
the density, especially, are all pure power laws and, for k < 1, they all diverge as x — > 0. As 
was mentioned in § I2.3[ this is not a serious problem because the singularity at the center 
can be removed by assuming that the central region rotates uniformly and that the density 
profile switches to a power law beyond a "core" radius x = x\. The core radius x\ can 
be chosen freely even for steep density profiles with k < — 1, but then the constant B in 
eq. ((27j) must be positive and it should be adjusted accordingly so that f(x) is everywhere 
positive and monotonically decreasing with x. It turns out that the physical requirement 



3 In the polytropic case, however, the analogue of eq. (|22|) . namely 

d f _ l dP \ 
dR \ P dhHl ) ~ ° ' 

is not satisfied by pure power-law profiles (except in the trivial case with P = constant); then, pressure- 
gradient variations do affect the structure of the underlying equilibrium solutions, but not in a dramatic 
fashion. 
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that f(x) > is weaker than the monotonicity condition that df/dx < 0. With A > to 
ensure that r(x) > and assuming that B > and k < — 1, we find from eq. (1271) that 
f{x) > if 

B ~ M • m 
where I = |fc + 1| > 0; and the stronger condition that df /dx < if 

^ > 4(l + f| • (31) 



ix\ 

Therefore, all composite models with a uniformly-rotating core region x < X\ must satisfy 
the stronger condition ( 13T1) for any choice of A > and /c < — 1 in their equilibrium density 
profiles (eq. [26]). 

The uniformly-rotating cores of the composite models discussed above call attention to 
another interesting feature: The density in these models must be constant and equal to (3$ 
in order to support this type of rotation (see also item 3 in the list of § !2.3p . So these models 
cannot obey the boundary condition that r(0) = 1 for centrally condensed structures. 
More generally, all the power-law solutions that we derived for the density in § 12.31 and 
all the composite models, although they are exact intrinsic solutions of the Lane-Emden 
equation with rotation (eq. [2U]), they do not solve the associated boundary-value problem. 
The question then is: How are the actual density profiles of centrally condensed equilibrium 
models going to behave once the proper set of boundary conditions {t(0) = 1, ^(0) = 0} 
are imposed at the center? 

The answer can be obtained by numerical integrations that enforce the desired central 
boundary conditions in eq. (120]) and in the analogous polytropic Lane-Emden equation with 
rotation; and by examining the analytic solution to the full boundary-value problem of the 
linear polytropic case with index n = 1 and uniform rotation. All the numerical (isothermal 
and polytropic) solutions, as well as the n = 1 analytic solution^ demonstrate routinely (see 
also Schneider & Schmitz [1995]) that the equilibrium density profiles lose their monotonicity 



4 For uniformly-rotating polytropic cylinders with n = 1, the solution to the boundary-value problem 
{t(0) = 1, £(0) = 0} is analytic (Robe 1968): 

t(x) = (1 - Pi) J (x) + 0l , 

where Jq(x) is the zeroth-order Bessel function of the first kind and all the other symbols are defined as in 
this work. (However, for n = 1 polytropes, c§ = 2Kp , where K = P/p 2 is the polytropic constant, and 
i?Q = K/2ttG.) It is clear in this solution that the Bessel function oscillates permanently about the r = 01 
line, which is a particular solution of the n = 1 linear ODE. 
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the first time they cross below the corresponding particular solutions (solutions analogous 
to those discussed in § 12.31 for the isothermal models). Once the first such crossing occurs 
at some relatively small radius (see Figs. [2] and [3] below), the actual physical solutions 
recognize the existence of the corresponding intrinsic solutions and they turn and oscillate 
permanently about the density level defined by these particular solutions. This of course 
happens because the particular solutions are fundamental "baseline" solutions of the ODE 
itself, regardless of externally-imposed conditions. When a set of external conditions are 
imposed at x = for physical reasons, the actual solutions emerging from the center do not 
match the baseline solutions (since the power-law behavior of the baseline is incompatible 
with the imposed conditions), but they are nonetheless attracted to them because the baseline 
solutions satisfy the ODE inherently. The result then is a permanent mismatch around the 
baseline that extends over all radii. This behavior is demonstrated in § 13.21 below, in the 
example model shown in Fig. [2] and, notably, in the composite isothermal model of the solar 
nebula shown in Fig. [31 

3. Application to Our Solar System 

The oscillatory behavior of the density profiles discussed in § 12.41 finds a natural appli- 
cation to the structure of the midplane of the solar nebula. For this application, we need a 
composite model as a baseline because such models allow for equilibrium density power laws 
that can be arbitrarily steep. To this end, we formulate in § 13.11 a composite model using 
the exact analytic solutions determined in § 12.31 above. Then, in § 13.21 we solve numerically 
the corresponding boundary-value problem that exhibits the same rotation profile as the 
analytic model; and we obtain an oscillatory solution for the density profile subject to the 
applied physical boundary conditions. Finally, we proceed to fit the density peaks of this 
model to the observed planetary distances in our solar system and we conclude in § 13.31 by 
determining important physical parameters associated with the structure, the dynamics, and 
the stability of the solar nebula. 

3.1. Composite Equilibrium Model 

For our baseline equilibrium model of the midplane of the solar nebula, we adopt a 
composite analytic solution in which the isothermal density profile has the form of a truncated 
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power law: 

{1 , if x < X\ 

(xjx) 6 , if xi < x < x 2 , (32) 

(xi/x 2 ) , if x > x 2 

where X\ is the radius of the constant-density core region, x 2 is the truncation radius beyond 
which the density remains constant at a low value, and the power-law index 5 is defined by 

5 = 1 - k 2 . (33) 

The condition that 5^2 implies that k ^ — 1 and excludes the logarithmic rotation laws 
from consideration: for k ^ — 1, logarithms are not introduced in the general form of the 
rotation law (eq. [2~7| ) by the lower branch of eq. (I2"gj) . Moreover, we consider below an even 
more limited range of indices, namely 5 > 2, since we are primarily interested in steep 
density profiles with k < — 1. 

Introducing X\ and x 2 in the above profile is equivalent to specifying three different 
values for the constant A in eq. f[2"oT) of § 12.31 These values are chosen so that the density 
profile remains continuous as it switches from one branch to the next, namely 

if x <x\ 

if X\ < x < x 2 . (34) 
if x > x 2 

The rotation law can then be determined from eqs. (T2TI) and by finding the values of the 
constant B that also make this profile continuous at the junctions where x = x\ and x = x 2 , 
namely 

if x < x\ 

if xi < x < x 2 . (35) 
if x > x 2 

As x — > x% from the right, the value of B in the intermediate branch satisfies marginally 
the monotonicity requirement (eq. [31]) determined in § 12.41 above. 

In practice, it is easier to integrate the differential equation ( |25l) in each of the three 
regions of the model and use the integration constants along with eq. fl34l to ensure continuity 



A 



x 



i • 



(xi/x 2 y 



B 



{xi/x 2 ) -{x l /x 2 ) 



51 
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at the junctions x\ and x 2 . Using either method, we find that the rotation law has the form 
1 , if x < xi 



/(*) 



1 



5 (xi/x)^ — 2 (xi/x) 



(xi/x 2 ) 2 - (x 1 /x 2 ) S ] (x 2 /x) 2 + (x 1 /x 2 f 



6-2 



if xi < x < x 2 



if x > x 2 



(36) 



It is easy to show that this rotation law obeys the physical requirements that f(x) > and 
df jdx < at all radii for any choice of the parameter set {x\ > 0, x 2 > x±, 5 > 2}. Notice 
that, outside the core, the rotation profile is monotonically decreasing everywhere; and that 
it becomes asymptotically flat at very large radii: as x — ► oo, then 

• (37) 



Thus, in contrast to the pure power-law density profiles (item 5 in § 12 . 31) . this composite 
profile exhibits nearly uniform rotation at very large distances. This is because the density 
of the model is not allowed to decrease at large distances; instead, it is kept constant at the 
low level shown in eq. (1321) for x > x 2 . 



The differential-rotation function f(x) of eq. ( 136]) and the corresponding equilibrium 
density profile Tb ase (x) / (3$ of eq. (132]) are shown in Figure 1 for x\ = 100, x 2 = 500, and 
for various choices of the power-law index 5 > 2. 



3.2. Solutions of the Boundary— Value Problem 
and Parameter Optimization to Planetary Distances 

The above composite equilibrium model is characterized by four free parameters: the 
core radius Xi, the truncation radius x 2 , the rotation parameter /3 < 1, and the power-law 
index 5 > 2. The density profile (eq. [32]) of this baseline solution of the Lane-Emden ODE 
(eq. |20j ) is not capable of satisfying physical boundary conditions at the center and it serves 
only as a mean approximation to the density of the corresponding physical model. The 
general form of the rotation law of the baseline (eq. [SB]) can, however, be adopted for the 
differential rotation f(x) of the physical model as well. Then eq. (132]) provides a prescription 
for the RHS of the Lane-Emden equation (T20]) which can thus be written as 

X d d / \ / \ 

— — x— lnr + r = r, (x) . (38) 
xdx dx base y 
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This ODE can be integrated numerically subject to the physical boundary conditions that 



An example is shown in Figure 2 for an equilibrium model with x\ = 200, X2 = 1000, 
(3q = 0.2, and 5 = 3. Notice that, in the linear scale of Fig. 2a, the actual density peaks 
(along the solid line) are approximately equidistant in regions where the baseline density 
(dashed line) is uniform; while in the intermediate region, they spread farther apart as they 
are trying to keep up with the steeply declining baseline. In the logarithmic scale of Fig. 2b, 
the same peaks appear to come closer together in the two areas where the baseline is flat 
and become equidistant in the middle area along the gradient of the baseline. 

A far more interesting application is shown in Figure 3. The Lane-Emden equation (I38p . 
subject to the physical boundary conditions fl39l . has been integrated numerically for various 
choices of the four free parameters. The resulting equilibrium profiles have been optimized 
for the present solar system assuming that their density maxima correspond to the observed 
semimajor axes of the planetary orbits out to and including Pluto. Not counting the central 
peak at x = 0, the third density maximum is always scaled to a distance of 1 AU during this 
nonlinear unconstrained optimization. In the best-fit model shown in Figure [3j the third 
density peak occurs at x = 44.564, implying that the length scale of the solar nebula in its 
isothermal phase was quite small (Rq = 0.022440 AU; see also § 13.31 below). 

The density maxima c?j of the best-fit model have been converted to AU and are listed 
in Table 2 along with the observed orbital semimajor axes a, taken from Table 1 above. The 
new dwarf planet Eris was not used in the optimization, but it is also listed at the bottom 
of Table 2 for comparison purposes. In addition, Table 3 displays all the extrema in the 
best-fitted oscillatory density profile out to 106 AU. Table 3 may be useful to observers who 
are trying to locate more dwarf planets in the outer solar system and to theorists who intend 
to build more sophisticated models of planitesimal accumulation in the solar nebula. 

Looking at orbital distances interior to the orbit of Eris in Table 2, we see that there 
are two peaks in the model, du = 48.73 AU and rf 12 = 58.27 AU, in which large dwarf 
planets have not been discovered. These two orbital distances lie beyond the outer "edge" or 
"gap" of the classical Kuiper belt (47-48 AU; e.g., Delsanti & Jewitt 2006), in an area where 
the number of orbiting objects decreases dramatically (e.g., Morbidelli, Brown, & Levison 
2003). Despite that, both peaks are located to within ~6% from two large objects: du is 
near Makemake (2005 FY 9 ), the third largest classical-Kuiper-belt object after Pluto and 
Haumea; and du is near 2002 TC302, the second largest scattered-disk object after Eris. 




(39) 
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Table 2: 

Locations of Density Peaks 
in the Best-Fit Model 
of the solar nebula 



Index 


Planet 


Semimajor 


Peak 




i 


Name 


Axis 


Location 


Error 






<H (AU) 


di (AU) 


(%) 


1 


Mercury 


0.387 


0.362 


-6.5 


2 


Venus 


0.723 


0.705 


-2.5 


3 


Earth 


1 


1 




4 


Mars 


1.524 


1.588 


4.2 


5 


Ceres 


2.765 


2.686 


-2.9 


6 


Jupiter 


5.203 


4.930 


-5.2 


7 


Saturn 


9.537 


9.843 


3.2 


8 


Uranus 


19.19 


20.09 


4.7 


9 


Neptune 


30.07 


29.66 


-1.4 


10 


Pluto 


39.48 


39.19 


-0.7 


11 




a 


48.73 




12 




b 


58.27 




13 


Eris 


67.89 c 


67.80 


-0.1 



Notes: 



(a) (136472) 2005 FY 9 is at a = 45.66 AU (d n deviates by 6.7%). 

(b) (84522) 2002 TC 302 is at a = 55.02 AU (d 12 deviates by 5.9%). 

(c) From Brown, Trujillo, & Rabinowitz (2005). 



The relative errors for each individual density peak of the model are also listed in Table 2. 
We see that the largest relative error is —6.5% for the first peak that corresponds to the 
orbit of Mercury. This deviation of ~0.025 AU is still quite small by solar-system standards. 
The optimization algorithm has minimized the mean relative error a for the first 10 orbits 
listed in Table 2. This was defined as a "standard deviation" by using the square deviations 
of 9 density peaks and excludes the peak that corresponds to the Earth's orbit which has 
zero deviation because of our scaling assumption that d% = 1 AU. Thus: 



a = 



N-2 



N 

E 



(di - eg)' 
a? 



(40) 



where N — 10 and the % — 3 term does not contribute to the sum (d 3 — a 3 = 0). By 
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Table 3: 

Density Extrema 
in the Best-Fit Model 
of the solar nebula 



Index 


Location 


Minimum 




Location 


Maximum 


% 


of Minimum 


Density 




of Maximum 


Density 




drain (AU) 


7~ (drain) 




dmax (AU) 


7" \d<max) 


1 


0.17901 


7.414 x 10" 


-2 


0.36162 


2.809 x 10" 1 


2 


0.53124 


1.069 x 10" 


-1 


0.70474 


2.466 x 10 _1 


3 


0.90231 


1.141 x 10" 


-1 


1 


1.216 x lO^ 1 


4 


1.4282 


3.571 x 10" 


-2 


1.5875 


3.757 x 10~ 2 


5 


2.4182 


9.483 x 10" 


-3 


2.6858 


9.821 x 10~ 3 


u 


A AQ^ 


9 D9^ v 1 IV 


-3 


zL (490^ 


o v 1 n -3 

Z.UOU A 1U 


7 


9.7814 


3.165 x 10" 


-4 


9.8431 


3.165 x 10~ 4 


8 


15.268 


1.586 x 10" 


-4 


20.093 


2.835 x 10" 4 


9 


24.855 


1.705 x 10" 


-4 


29.661 


2.714 x 10^ 4 


10 


34.421 


1.773 x 10" 


-4 


39.193 


2.642 x lO" 4 


11 


43.962 


1.818 x 10" 


-4 


48.726 


2.593 x 10" 4 


12 


53.490 


1.851 x 10" 


-4 


58.266 


2.557 x 10~ 4 


13 


63.028 


1.876 x 10" 


-4 


67.800 


2.529 x 10" 4 


14 


72.557 


1.897 x 10" 


-4 


77.336 


2.507 x 10~ 4 


15 


82.092 


1.914 x 10" 


-4 


86.857 


2.488 x 10~ 4 


16 


91.626 


1.928 x 10" 


-4 


96.399 


2.473 x 10" 4 


17 


101.15 


1.940 x 10" 


-4 


105.91 


2.460 x 10~ 4 



construction, this definition of the mean relative error places more weight to the orbital 
distances of planets near the Sun and allows for larger errors in the locations of the outer 
density peaks. Despite this skewing of the fit, the mean relative error for the best-fit model 
is a = 4.1%, affirming that our simple equilibrium model succeeds in matching all of the 
observed planetary distances very well. 

3.3. Physical Parameters of the Solar Nebula 

The parameters of the best-fit model determined by the optimizing algorithm along 
with the scaling assumption that d 3 = 1 AU (see Table 2) constitute a set of important 
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dynamical parameters of the long-gone solar nebula: 

0o = 0.41465 
S = 2.5362 
k = 1 - 6 = -1.5362 



x 2 



R = 0.022440 AU 
= 36.511 := 0.81931 AU 
= 505.45 := 11.342 AU 



(41) 



The value of the rotation parameter (3q ~ 0.4 indicates that the rotation of the isothermal 
nebula was moderate, about 40% of the maximum value allowed by self-gravity. This result 
also indicates that the best-fit model (and the other models discussed here) is stable to 
nonaxisymmetric disturbances because (3q is much below the critical value of /3* = 0.7. 
This critical value can be obtained from the a-parameter criterion for stability of rotating, 
self-gravitating, gaseous systems (a < 0.35; Christodoulou, Shlosman, & Tohline 1995) by 
combining the definition f3 = Q Q /Qj with the definition of a for disks where a = VL$/2Vtj to 
get the relation /3 = 2a which implies that (3* = 0.7. Of course the models discussed in this 
section are stable to axisymmetric disturbances as well, since they all satisfy the Rayleigh 
criterion. 

The value of the power-law index 5 ~ 2.5 indicates that the density profile of the 
differentially-rotating region of the nebula declined with radius, on average, as R~ 2 5 . This 
region extended from a radius of x\ ~ 0.8 AU out to a radius of x 2 ~ 11.3 AU. Of these 
parameters, only x 2 is slightly uncertain because the mean relative error a (eg. [40]) places 
less weight to the orbits of the outer, nearly equidistant planets which, in effect, determine 
the truncation radius. As a result, the optimization procedure also finds additional "near- 
minima" of high quality (cr = 4.1%-4.7%), among which the most extreme model has the 
following parameter values: 5 = 2.5040, fa = 0.3806, #0 = 0.02032 AU, x x = 40 
(:= 0.81 AU), and x 2 = 575 (:= 11.7 AU). A comparison between these values and the 
best-fit values listed in eq. (j4"Tl) gives us an idea about how shallow the region around the 
true minimum is in the four-parameter space of the model. Note, in particular, that the 
power-law index does not differ from 5 = 2.5 by more than 1.5% and the rotation parameter 
does not differ from (3q = 0.4 by more than 5% in any of the high-quality fits to the data. 

The cylindrical Lane-Emden equation that we have solved can serve as a model of 
differentially-rotating disks supported by thermal pressure in the ^-direction, so we expect 
that the scale height from pressure support will be H oc R down the radial density gradientB 



5 The vertical scale height of a pressure-supported disk is H ~ cq/Q ~ Rcq/v, where v is the rotation 
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In this case, the value of 5 = 2.5 implies that the corresponding power-law index in the 
surface-density profile (S oc R~ s+1 ) of the nebular disk was k — 1 — S — —1.5. This value 
is virtually identical to that obtained by Weidenschilling (1977) who derived an estimate 
of the surface density distribution of the protoplanetary disk by smearing out the observed 
planetary masses over annular rings in the disk's midplane and then applied a correction 
to this distribution by adding the appropriate amount of volatiles in order to bring the 
elemental abundance of the gas up to solar composition. Our analytical work and modeling 
effort also provide a direct approach to the same problem, but from a different angle than 
that conceived by Weidenschilling. The unambiguous congruence of the results obtained by 
these two disparate methods is rather astonishing and suggests strongly that the surface 
density profile T,(R) of the solar nebula was indeed exhibiting an R~ 15 power law in the 
isothermal phase of its evolution. 

Furthermore, our work also helps in delineating the fundamental physics behind such 
a mean surface density profile in the midplane of the solar nebula. With the aid of our 
best-fit model, we can deduce substantial new information concerning the structure and the 
dynamics of the nebular disk. In addition to the structural and rotation parameters discussed 
above, we can use our analytic baseline model in order to probe the dynamical state of the 
protoplanetary disk in its isothermal phase as follows. 



3.3.1. Equation of State 

Using the length scale of the disk (R = 0.022440 AU) in eq. ([9]), we can write an 
equation of state of the form 

c 2 

= iwGR 2 = 9.45 x 10 16 cm 5 g" 1 s" 2 , (42) 

Po 

where Co and po m ay be thought of as the local sound speed and the local density in the 
inner disk, respectively. For an isothermal gas at temperature T, Cg = 1ZT/JI, where Jl 
is the mean molecular weight and 1Z is the universal gas constant. Hence, eq. (14*21) can be 
written in the form 

po = 8.80 x 10" 10 gem" 3 , (43) 



speed. In our models, v is asymptotically flat when the density exhibits a power-law profile and cq is 
constant, leading thus to the approximate relation H oc R. On the other hand, the scale height is 
H ~ co/ilo = V^Ro/Po = 3.4i?o in the core of our best-fit model. As a result, the mass estimates 
calculated in § 13.3.31 for a disk with H = R are too conservative. 
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where T and p are measured in degrees Kelvin and g mol , respectively. 

For a cold disk of gas with T = 10 K and p = 2.34 g mol -1 (molecular hydrogen and 
neutral helium with fractional abundances X = 0.70 and Y = 0.28 by mass), we find that 

p = 3.76 x 10" 9 g cm" 3 . (44) 

This value is comfortably larger than the well-known threshold for planet formation in the 
solar nebula (p* ~ 10~ 9 g cm -3 ; see e.g., Lissauer [1993]) and implies that the conditions for 
planet formation were already in place, at least in the inner disk, in the isothermal phase of 
the disk's evolution. 



3.3.2. Rotational State 

Using the characteristic density of the inner disk (eq. [44] ) in the definition of Qj = 
y/2iiGpo, we can determine the Jeans frequency of the disk: 

ttj = 3.97 x 10~ 8 rad s" 1 . (45) 

Then, using the value (3q = 0.41465 (eq.|41j) in the definition of (3q = flo/flj, we can 
determine the angular velocity of the uniformly-rotating core (i?i < 0.81931 AU): 

tt = 1.65 x 10~ 8 rad s" 1 . (46) 

For reference, this value of Qq corresponds to an orbital period of 12 yr. In the present solar 
system, that would correspond to a Keplerian orbit with semimajor axis a = 5.24 AU. So 
the core of the nebula was rotating about as slowly as Jupiter is currently revolving around 
the Sun. 

Finally, the constant asymptotic value of the angular velocity for x >> x 2 is 

also a characteristic rotational parameter of the nebula (because the outer region of uniform 
density is necessary for the formation of the nearly equidistant outer planets). Using eqs. (j3j) 
and (|37|) . we find that 

= Vl (^-j = 5.89 x 10~ 10 rad s" 1 , (47) 

where the values listed in eq. ( |4"T1) and eq. (1461) were used in the numerical evaluation. This 
value of Qoo is 28 times smaller than Qq; and corresponds to an orbital period of 338 yr 
and to a Keplerian orbit with a = 48.5 AU near the outer edge of the classical Kuiper belt 
in the present solar system. 
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3.3.3. Mass Distribution 

Assuming that our model disk is uniform in the z-direction over a thin layer of thickness 
2Ro, we can multiply its mass per unit length by this factor and we can define the total 
mass near the midplane of the solar nebula out to a maximum radius R max as 

M = 2R ■ [ 2np(R)RdR = 2M • / r{x)xdx , (48) 
Jo Jo 

where the constant Mq is given by 

M = 27rp #o = 4 - 47 x 10_7 M e , ( 49 ) 

where eq. ( 1441) and Rq = 0.022440 AU were used in the numerical evaluation. M is 
approximately the mass contained to within one length scale from the center of the disk and 
implies that the central surface density is 

S = = 2R -p « 2520 gem' 2 . (50) 

Using this value, the core radius Ri = 0.81931 AU, and the power-law index k ~ —1.5, we 
can write the surface density profile of the solar nebula for R\ < R < R2 in the form 

E(R) = 1870 (rjw) g cm " 2 ' ( 51 ) 

where R is measured in AU. The value at 1 AU is lower than that estimated by Weiden- 
schilling (1977); but it agrees very well with Hayashi's (1981) competing result that was 
obtained by the same method and with the same data, and also with Kuchner's (2004) 
average estimate obtained from a similar analysis of 26 multiple-planet extrasolar systems. 

Adopting now the baseline density profile (eq. [32]) as an approximation to the actual 
density distribution and using the parameters listed in eq. (jHJ), we can evaluate the integral 
of eq. (j!5|) over the three regions of the baseline out to, e.g., R m ax = 50 AU (x max = 2228): 



M = 2M Q fil 



asi rx 2 / Xx \ft f Xmax f X\\ 

xdx + / ( — J xdx + I — ) xdx 

J xi ^ x ' Jx 2 \ x 2/ 

1 x 10 -4 M Q + 3 x 10" 4 M & + 5 x 10 -4 M (52) 
0.001 M . 



The total mass is one order of magnitude smaller than the low end of the estimate for the 
" minimum-mass solar nebula" (0.01 M Q ; Weidenschilling 1977). This is not surprising since 
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we have adopted a very small value (Rq) for the vertical scale of the disk. In the classical 
scenario of cloud collapse, this mass near the midplane will be enhanced substantially as 
more matter from the parent cloud will continue falling onto the disk (see § 13.3.41 below for 
more details). It is however interesting that the inner disk already has high enough densities 
(eq. [44]) to begin the process of planet formation so early in its evolution, before the central 
protostar grows to become gravitationally dominant or radial accretion becomes important 
in the gas. 



3.3.4- Integral Properties of the Core 

Integrating over the mass distribution of the core of the baseline model, we can esti- 
mate important dynamical properties, such as the core mass and angular momentum. The 
analytic estimates can then be compared to the corresponding results from protostellar col- 
lapse simulations. Such simulations have been recently reviewed by Tohline (2002) who also 
summarized the main results over which there seems to be widespread agreement among 
researchers working in the field for almost 40 years (see the discussion of Tohline centered 
around his Figure 2 and Table 2). 

The mass of the core found in § 13.3.31 above. M± ~ 10~ 4 M , is 40 times smaller than 
the typical mass of a collapsing cloud core in which the Jeans instability may be temporarily 
halted by thermal pressure and an adiabatic quasi-equilibrium may then be established (see 
case B in Table 2 of Tohline [2002]). This means that the core of our quasi-equilibrium 
disk will remain isothermal and it will continue to accumulate mass from the infalling cloud 
beyond the point described here. 



Furthermore, the total angular momentum of the core of the model is 

-Hi 

'o Jo 
This implies that the specific angular momentum of the matter in the core is 



pRi r%i i 

Li = 2 Rq ■ / n R 2 ■ 2irp RdR = 2M n R 2 ■ / x 3 dx = -M tt R 2 xf . (53) 
Jo Jo 2 



Mi 1 1 2(3 2 



7.19 x 10 18 cm 2 s _1 , (54) 



which is about one-half of the corresponding estimate given by Tohline (2002) for a cloud 
core at the endpoint of its isothermal evolution. Clearly, in a realistic setting, this type of 
low-mass core has the potential to grow by accreting a lot more matter of low specific angular 
momentum, the kind that can settle on to the central region from the vertical direction. 
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4. Summary and Discussion 
4.1. Summary 

We have presented exact analytic and numerical solutions of the axisymmetric, cylin- 
drically symmetric, isothermal Lane-Emden equation with and without rotation (eqs. [20] 
and [10], respectively). This second-order ODE describes the radial equilibrium configura- 
tions that are available to a self-gravitating perfect fluid in which the vertical variation of 
the enthalpy gradient, namely 

d_ ( _idP\ 
dz \ dz J ' 

is negligible. The enthalpy gradient p~ 1 dP/dz, where p is the density and P is the pressure, 
is the term that establishes vertical hydrostatic equilibrium when it successfully balances the 
gravitational force d^/dz in the z-direction, but it is its vertical variation that is actually 
ignored when the cylindrical Lane-Emden equation is considered. This term is zero identi- 
cally for a cylindrical model of infinite vertical extent; such a model may be applicable to 
elongated, filamentary, star-forming regions in which z > R. In other astrophysical applica- 
tions, especially those dealing with gaseous disks that concern us in this work, this second 
derivative is not identically zero everywhere; but it is ignored on the basis that it vanishes 
on the symmetry plane and hopefully it remains negligible away from that plane over some 
layer of the astrophysical disk (as in the exact models of Schmitz & Ebert [1986]). In the 
specific case of the solar nebula, this assumption is probably valid because the disk is not in 
vertical hydrostatic equilibrium and vertical pressure support is smaller than the weight of 
the infalling matter everywhere, except near the midplane of the disk where the gradients 
tend to zero by symmetry anyway. 

The rotating analytic solutions are exact intrinsic solutions of the Lane-Emden equa- 
tion (120j) . Their density profiles are all pure power laws or combinations of power laws 
(eqs. [26] and |32j). These solutions are determined solely by the balance between the gradi- 
ents of the centrifugal and the gravitational accelerations (eq. [29]) because for an isothermal 
fluid with a power-law density profile, the radial variation of the enthalpy gradient, namely 

d ( _ x dP \ 
dR V dlnR ) ' 

is zero identically (see eq. [22] and the discussion at the end of § 12. 3p . The differential- 
rotation profiles that can support such power-law densities are quite general and varied 
(eqs. [27] and [5H]). They contain three integration constants that can be used to compose 
models of equilibrium disks in which the density profiles are nonsingular and arbitrarily steep 
in radius. We have created one such composite equilibrium model and we have applied it to 
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the structure of the midplane of the solar nebula (§ [3]). In this model, we have connected 
a central core region of constant density and an outer region of constant density with a 
power-law profile. The rotation parameter and the size of the core as well as the power-law 
index and the truncation radius of the intermediate region are the four model parameters to 
be determined based on the presently observed solar system. Beyond the truncation radius, 
the density remains constant at a low level. This turns out to be the necessary nebular 
background for the formation of nearly equidistant outer planets and its physical origin 
deserves further investigation. 

The density profile of the composite equilibrium disk model is no longer an exact solution 
when physical boundary conditions are imposed at the center of the disk. The exact solution 
to the associated boundary-value problem is obtained by numerical integration that enforces 
the proper central boundary conditions in a model with the same differential-rotation law 
as the composite analytic model. We have found that the intrinsic analytic density profile is 
a good approximation to the actual numerical solution, i.e., it is a "baseline" that exhibits 
on average all the important features of the exact solution to the boundary-value problem, 
except one: the actual density profile is permanently oscillatory in radius. In fact, it oscillates 
around the baseline solution because it is attracted to this inherent particular solution of the 
Lane-Emden equation, but the two solutions cannot match in any finite segment since such 
coincidence is strictly prohibited by the imposed central boundary conditions (see Fig. 2 and 
the discussion in § 12.41) . 

The radial density peaks of the actual solution to the boundary- value problem corre- 
spond to local minima of the gravitational potential in the midplane of the protoplanetary 
disk. In a cloud collapse scenario, the density of the gas and the concentration of conden- 
sates (dust, rocks, and ices) will be enhanced inside these local gravitational potential wells 
as more matter rains down onto the disk from the surrounding protostellar cloud. It is not 
unreasonable to consider that the added material will, in turn, deepen further these potential 
wells which may thus become feasible sites of planet formation. Based on this picture, we 
have proceeded to fit the density peaks of the exact numerical solutions of the Lane-Emden 
boundary-value problem to the presently observed locations of the planets around the Sun 
and the satellites of the Jovian planets. Our best-fit model for the planets in our solar 
system (Fig. 3 and Tables 2 and 3) is described in § 13.21 and the dynamical parameters 
obtained for the solar nebula are described in § 13.31 The mean relative error for this best-fit 
model is 4.1% and the relative deviations for individual planets do not exceed 6.5%. Fur- 
thermore, this model suggests that the radial surface density profile of the protoplanetary 
disk between roughly 0.8 AU and 11.3 AU declined as R~ 15 , a result that agrees fully with 
the classic determinations of Weidenschilling (1977) and Hayashi (1981) which were based 
on the presently observed planetary masses. We find the confluence of these very diverse 
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approaches encouraging and supportive of the view that the fundamental notions behind 
our model are indeed relevant to the problem of the radial distribution of planetary orbits. 
We would also like to note that the same type of modeling works very well for the regular 
satellites of Jupiter and for the 5 planets of 55 Cancri, and we plan to present these results 
in future publications. 



4.2. Discussion 



Our work is the first to produce exact solutions of the nonlinear isothermal Lane- 
Emden equation with differential rotation that exhibit a pronounced oscillatory behavior, a 
feature attributed exclusively to the proper choice of the central boundary conditions. The 
roughly arithmetic or geometric spacings of the corresponding density maxima depend on 
the local slope of the differential rotation profile and provide a transparent, physically-based, 
and soundly formulated reason for producing a sequence of distinct density enhancements 
at radial positions that agree with the present planetary distances, given that the latter 
are observed to follow an arithmetic, then a geometric, and finally again an arithmetic 
progression (see Table 1). Hence, it should be remarked that we now have at hand, for the 
first time, a plain explanation of the so-called Titius-Bode "law" of planetary distances. 
This empirical algorithm has been known for 240 years but with no underlying physical 
justification (see e.g. Graner & Dubrulle [1994] and Hayes & Tremaine [1998]). Our best- 
fit model of the midplane of the solar nebula (Fig. 3 and Table 2) shows a progression 
of density enhancements that are stretched farther apart from one another, as the actual 
density profile decreases sharply in trying to keep up with the steeply declining baseline. The 
relative spacing of the peaks is roughly geometric between 0.7 and 20 AU, and this explains 
the success of the Titius-Bode algorithm for a large number of planetary orbits. The change 
of the underlying density profile to constant beyond a certain distance in our model can then 
also explain why this algorithm fails at large distances. A similar account of the nongeometric 
spacing of Mercury also finds a straightforward explanation in terms of a uniform central core 
in our model. So, just as many researchers have suspected in the past (see § 11.21) . the Titius- 
Bode algorithm has no underlying physical principle underneath its phenomenology; the 
algorithm has been successful only to the extent that it has managed to exploit empirically 
the most pronounced feature of the radial density profile of the solar nebula (the roughly 
geometric spacing of the peaks between 0.7 and 20 AU). Finally, we should point out that 
the above resolution of the long-standing problem of planetary "order" in our solar system 
is very anticlimactic, as the long-sought physical explanation turns out to lie entirely within 
the realm of conventional physics and does not need to invoke unrestrained numerology of the 
Titius-Bode type, exotic physics, "new" dynamical laws, or arbitrary "universal" constants 
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and ad-hoc solar-system "quantizations" such as those proposed in the past. 

We view this work only as a first (but compelling) step toward understanding the sys- 
tematic formation of planets; the evolution that results in the formation of the entire solar 
system likely involves additional processes and physics (e.g., dissipative disk accretion, mag- 
netic fields) that are necessary in delineating its present structure and composition. However, 
we believe that the condensed solid cores inside the gravitational potential minima of the 
solar nebula could possibly survive the subsequent evolution and could form planets at the 
same orbital locations. Furthermore, we hope that the results presented here will provide 
the motivation behind the formulation of some novel hypotheses concerning the formation 
and the dynamical evolution of certain sectors in the early solar system and in the recently 
discovered extrasolar planetary systems (see the reviews of Lissauer [1993], Beckwith Sz Sar- 
gent [1996], and Marcy et al. [2005]). Some of the ideas that emerge from our study are 
outlined below. 



4-2.1. The Distribution of Condensates in the Protoplanetary Disk 

Our results support the segregation of condensates at specific radii in which the gas 
density is larger than the mean background value (Fig. 3 and Table 2). These orbital 
locations in the midplane of the disk are ideal sites for growing planitesimals by accumulation 
of smaller bodies in a systematic (nonrandom) way. In the best -fit density profile r oc x 
with 5 = 2.5, the absolute spacing of these sites of concentrated material depends on a single 
parameter, the length scale of the disk (Rq = 0.022440 AU) which is a measure of the entropy 
content of the gas (see eq. [4"2"]). 

The segregation of condensates in initially shallow (see Table 3) gravitational poten- 
tial wells should be tested by computational experiments because this kind of work may 
potentially lead to improved protoplanetary models. Numerical work is quite common in 
investigations of the early solar system but the simulations have always started from arbi- 
trary initial conditions (e.g., Wetherill 1989; Ruzmaikina, Safronov, & Weidenschilling 1989; 
Bodenheimer, Ruzmaikina, & Mathieu 1993; Boss 1995) because of the limited information 
that can be collected from observations, meteorites, and laboratory experiments. Our results 
can potentially provide a better handle on the appropriate initial conditions for simulations 
of this type. 
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4-2.2. Planet Migration in Our Solar System 

Numerical investigations of planet migration have intensified over the past 10 years 
(see the progress report by Levison et al. [2007]) in an attempt to understand the complex 
dynamics observed in the Kuiper belt: the trapping of some small objects in the mean- 
motion resonances of Neptune, the gravitational scattering of other small objects to high 
inclinations, the severe mass deficiency (by at least a factor of 100) of the entire Kuiper 
belt, and the abrupt outer "edge" or "gap" of the classical Kuiper belt at 47-48 AU (e.g., 
Delsanti & Jewitt 2006). Our results do not support the idea that the gaseous giant planets 
have secularly migrated over time to the presently observed radial locations (Fernandez & 
Ip 1984; Malhotra 1993, 1995; Gomes, Morbidelli, & Levison 2004; Gomes et al. 2005). 
The quality of our best-fit model of planetary distances is so high (see § 13 . 2j) that it seems 
improbable that any planet that scattered small planitesimals has managed to jump out of 
its local gravitational potential well and has moved to a new location inside an adjacent 
potential well. 

On the other hand, the displacement of a planet within its potential well appears viable 
in the context of our model, and the potential wells of the gaseous giants are quite large in 
size (see Table 3 and Fig. 3a): Uranus and Neptune are both orbiting on potential minima of 
radial half-width Ad ~ ±5 AU; Saturn's minimum has Ad ~ +5/ — 0.06 AU; and Jupiter's 
minimum has Ad ~ +5/ — 0.4 AU. These values are not very different than the migration 
distances used so far in the numerical simulations; only Saturn presents a challenge because 
its potential well is too steep at smaller radii to allow for an outward migration of more than 
~0.1 AU to its present location. Perhaps new models of planet migration can be constructed 
in which these results will be taken into account. 



4-2.3. An Inner Gap in Protostellar Disks 

Fig. 3b and Table 3 show that the first density minimum in the best-fit model of the solar 
nebula occurs at d = 0.179 AU. This suggests that there were no likely planetary sites interior 
to the orbit of Mercury. It also suggests that there was a significantly lower concentration of 
condensates within this minimum, in the area below the baseline that extends approximately 
from 0.1 AU to 0.3 AU. Since we expect that our solar system is in no way special but rather 
representative of large planetary systems around solar-type stars, we believe that the same 
deficiency of solids should also exist in the inner regions of other protostellar disks that are 
currently in the process of forming their protostars and protoplanets. We think that such an 
inner gap in condensates is currently inferred in the circumstellar disks of some pre-main- 
sequence stars; these disks show almost no near-infrared radiation emanating from the inner 
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0.2 - 0.3 AU (Strom, Edwards, & Skrutskie 1993; Beckwith & Sargent 1993; Millan-Gabet 
et al. 2001; Akeson et al. 2005; Fedele et al. 2008). The alternative view is that an 
orbiting protoplanet may have cleared a gap in that area. Although this may be possible 
in small heavy disks where large planets could form very close to their stars (Christodoulou 
& Kazanas 2008), it cannot explain large systems like our own in which gas giants are not 
expected to form within the inner 0.3 AU. 

4-2.4- Extrasolar Planetary Systems 

It is believed on theoretical grounds that terrestrial planets are commonly formed near 
their central stars and that it is unlikely to have a gaseous giant form closer to the center than 
about 4 AU (Boss 1995). Accordingly, our solar system is presumed to be typical of planetary 
systems around other stars. On the other hand, most of the extrasolar planets found to 
date do not fit in this theoretical picture. The currently available sample of extrasolar 
planets is strongly biased toward giant planets because small perturbing masses cannot be 
easily resolved by observations, so it is understood that it will take more time before large 
extrasolar systems comparable to our own can be detected. Nevertheless, there is presently 
no doubt that in many extrasolar systems, a large number of massive planets exist very close 
to their stars (the trail of discovery started with the planet at 0.05 AU in 51 Peg found 
by Mayor & Queloz [1995] and the planet at 0.48 AU in 70 Vir found by Marcy & Butler 
[1996]). This discrepancy seems to point to the hypothesis that the cores of these planets 
have formed farther away and that the objects have migrated inward to the observationally 
inferred distances (Kary & Lissauer 1995; Lin, Bodenheimer, & Richardson 1996), perhaps 
destroying in the process the orbits of terrestrial planets that could have existed closer to 
central stars. 

Migrating giant planets cannot be accounted for in our model, so we will have to wait 
until observations can detect some large extrasolar systems similar to our own. However, 
two small multiple-planet systems have already been detected which are different than our 
solar system, but their planetary distributions appear to be more in line with a sequence of 
well-ordered density peaks as predicted by our model. These systems are 55 Cancri (Fischer 
et al. 2008) and HD 37124 (Vogt et al. 2005). Using updated data from the "Catalog 
of Nearby Exoplanets" ( |http: //exoplanets.org| and Butler et al. [2006]), we find that the 
ratios of semimajor axes of the three innermost planets in these two systems are 02/01 ~ 3 
and 03/02 — 2. These ratios are barely larger than those in our solar system (1.9 and 
1.4, respectively). In contrast, all the other systems in the Catalog with three or more 
exoplanets have 02/01 > 6 and/or 03/02 > 3, ratios that are too large compared to those 



-35 - 



in 55 Cancri, HD 37124, and our solar system. 

We see then that the relative distributions of planetary orbits in 55 Cancri and HD 37124 
are effectively the same, despite the fact that these two well-ordered systems are much smaller 
than our solar system and their central stars have widely different chemical compositions 
(55 Cancri is very metal-rich, while HD 37124 is very metal-poor). 55 Cancri is now believed 
to have 5 planets in remarkably circular orbits (Fischer et al. 2008), and this number is 
enough to allow for our type of modeling for this system. We have carried out the analysis 
of 55 Cancri and we present our results in a companion paper (Christodoulou & Kazanas 
2008). 
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FIGURE CAPTIONS 

Fig. 1. — Analytic density and rotation profiles of the composite equilibrium models de- 
scribed by eqs. (152]) and (!36l) for xi = 100, x 2 = 500, and 5 = 2.5, 3, and 4. The density 
profile Tb ase (x) / (3$ is uniform for x < X\ and for x > x 2 ; and it follows the power law x~^ 
in the in-between region. The rotation profile f(x) is uniform for x < x± and monotoni- 
cally decreasing for x > x±, as specified by eq. ( 136]) : at very large radii (for x » X2), f(x) 
approaches the constant asymptotic value 

Fig. 2. — Equilibrium density profile for a model with rotation parameter /5 = 0.2 and a 
composite rotation profile given by eq. (136]) with x\ = 200, x 2 = 1000, and 5 = 3. Frame 
(a) shows the radial distance x on a linear scale out to 1 = 2000. Frame (b) shows the 
same radial distance on a logarithmic scale out to In a; = 8 (x = 2981). The physical 
density r(x) (solid line) satisfies the boundary conditions (13"9~|) and, as a result, it is forced to 
oscillate permanently about the inherent baseline solution r base (x) (eq. [52] , dashed line) of 
the Lane-Emden ODE fl38l) . The nonrotating analytic solution (eq. [17] . dash-dotted line) 
is also shown for reference. 

Fig. 3. — Equilibrium density profile for the midplane of the solar nebula. The composite 
model described in § 13.11 (eq. [32], dashed line) has been adopted for the RHS of the Lane- 
Emden equation (1381 and this ODE has been integrated numerically subject to the physical 
boundary conditions (139]) . The resulting solution (solid line) has been fitted to the present 
solar system so that its density maxima (dots) correspond to the observed semimajor axes 
of the planetary orbits (open circles). The third density maximum is always scaled to a 
distance of 1 AU in this procedure; in the best-fit model shown here, it occurs at x = 44.564 
and implies that the length scale of the solar nebula was Rq = 0.022440 AU. The mean 
relative error of the fit is 4.1%, affirming that this simple equilibrium model produces an 
incomparable match to the observed data. Frame (a) shows the radial distance d on a linear 
scale out to d = 40 AU. Frame (b) shows the same radial distance on a logarithmic scale 
out to \nd = 4.5 (d = 90 AU). The nonrotating analytic solution (eq. [T7] . dash-dotted 
line) is also shown for reference. 
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